// Copyright (C) 2009,2010,2011,2012  Marco Restelli
//
// This file is part of:
//   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
//
// LDGH is free software: you can redistribute it and/or modify it
// under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// LDGH is distributed in the hope that it will be useful, but WITHOUT
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
// or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
// License for more details.
//
// You should have received a copy of the GNU General Public License
// along with LDGH. If not, see <http://www.gnu.org/licenses/>.
//
// author: Marco Restelli                   <marco.restelli@gmail.com>

// Set the the domain size
Lx = 220.0e3;
Ly =  30.0e3;
h  =   2.0e3;
a  =  10.0e3;
xm = 10.0*a;

// Set the the mesh parameters
NPx = 161;
NPy = 113;
dx = Lx/NPx;
x_ref = 125.0e3;
y_ref = 12.5e3;
sigmax = 40;
sigmay = 40;
Dmin = 0.1*dx;
regular = 1;

// Build the mesh
Nx = 3*NPx; // make sure the mountain is well represented
For t In {1:Nx}
  x = Lx/(Nx-1)*(t-1);
  x0 = x - xm;
  y = (h*a^2)/(x0^2+a^2);
  Point(t) = {x,y,0.0,dx};
EndFor
Spline(1) = {1:Nx};

Point(Nx+1) = { Lx,Ly,0.0,dx};
Point(Nx+2) = {0.0,Ly,0.0,dx};

Line(2) = {  Nx,Nx+1};
Line(3) = {Nx+1,Nx+2};
Line(4) = {Nx+2,   1};

Line Loop(5) = {1,2,3,4};
Plane Surface(6) = {5};

If(regular == 1)
  Transfinite Line{1,3} = NPx;
  Transfinite Line{2,4} = NPy;
  Transfinite Surface{6} = {1,Nx,Nx+1,Nx+2};
EndIf

// Localized refinement
Field[1] = MathEval;
Field[1].F = Sprintf(
  "sqrt( ((x-%g)/%g)^2 + ((y-%g)/%g)^2)+%g",
  x_ref , sigmax , y_ref , sigmay , Dmin );
Background Field = 1;

